Many-body resonant tunneling in the Wannier-Stark system 
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We study an experimentally realizable paradigm of a complex many-body quantum system, a two- 
band Wannier-Stark model, for which diffusion in Hilbert space as well as many-body Landau-Zener 
processes can be engineered. A clean cross-over between regular to quantum chaotic spectra is found 
within the many-level avoided crossings at resonant tunneling conditions. The spectral properties 
determine the evolution of states across a cascade of Landau-Zener events, which can be used to 
study diffusion and to dynamically prepare target states interesting for experimental applications. 
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Introduction. The rapid development of high precision 
techniques for the experimental control of ultracold quan- 
tum gases offers a clean way to study static and dynam- 
ical properties of interacting many-body lattice systems 
[TJ [2] . Of particular interest are realizations of strongly 
correlated or complex systems of many particles. As- 
pects of integrability of such quantum systems are cru- 
cial for predictions, for instance, of their thermalization 
[3JH]. Moreover, experiments controlling the populations 
of higher orbitals and bands in periodic potentials are at 
hand now [SJ [6] , which allow for the study of complex 
many degree-of- freedom systems. 

In this Letter we present a case study of a complex 
many-body system including two strongly coupled en- 
ergy bands. As we sketch below it can readily be re- 
alized experimentally with ultracold bosons in acceler- 
ated (or equivalently tilted) optical lattices [7 . Our 
system is a generic model for which the dynamics can 
be steered by the parameters (tilt, interaction strength, 
potential depth), thus implementing very different dy- 
namics. While the mean- field transport (weakly interact- 
ing limit) is well known for this system [7J [8], even the 
simplest many body version, a one-band Bose-Hubbard 
model with(out) tilt, allows one to change between reg- 
ular and quantum chaotic evolutions [9HT21. Our system 
is a paradigm for which, depending on the choice of pa- 
rameters, full complexity in the interband transport and 
thermalization in Hilbert space can be reached. 

Strongest interband coupling occurs at resonantly en- 
hanced tunneling (RET) between the two bands [SJ [13] , 
where we find a crossover from regular to quantum 
chaotic spectral statistics with increasing density of 
states. These complex energy spectra determine the 
transport across the many avoided crossings at RET 
when the force F is varied. This not only generalizes 
results on the weakly interacting limit [14], but relates 
to the largely open problem of many-body Landau-Zener 
processes in the precence of strong particle interactions 
[TTJ [15] [16]. As two direct applications we show how 
the spectral properties influence the diffusion in Hilbert 
space, and, moreover, how to prepare target states with 
a high fidelity in the upper band by controlled sweeps of 
F(t) through the RET regime. 




FIG. 1: (Color online) The system: (a) Many-body processes 
of the two-band Bose-Hubbard Hamiltonian for a bichromatic 
one dimensional tilted optical lattice, (b) Resonant Enhanced 
Tunneling (RET) condition for the nearest neighboring dou- 
ble wells, i.e, for a first order resonance. 



The two-band system. Single-particle energy 
bands are nowadays routinely engineered in ex- 
periments using a bichromatic optical lattice 
V(x) = — Vo [cos(2/c^x) + sq cos(4/c^x + </>)], being 
the recoil momentum [17 . Choosing, e.g., the relative 
phase (j) = 7r and appropriate values of sq — Vd/Vo 
we have a situation for which the lowest two energy 
bands are well separated from all higher bands [18 . An 
additional Stark force stimulates the quantum transport 
along the lattice [7J [8] [10] [TTJ [14] and, at the same time, 
it couples the bands. In the single particle picture, this 
coupling is maximal at specific tilts F r w A^/27rr, 
being the energy gap between the two bands. At those 
values RET occurs between levels distancing r double 
wells, see Fig. [TJ and the integer r is called the order 
of the resonance [8]. The corresponding many particle 
problem can be described in the tight-binding limit by a 
two-band Bose-Hubbard Hamiltonian 
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FIG. 2: (Color online) The spectrum: (a) Wannier-Stark lad- 
ders from two bands for the single particle case, with gap 
A g « 1.61 at Vb = 5 and so = 2. Avoided crossings appear 
at the resonances F r =i = 0.1934 with width A r <C A g . (b) 
The many-body spectrum revealing the local complex struc- 
ture of many avoided crossings even for this relatively small 
system (N/L = 5/3). £o(F) is the F dependent center of the 
FZ. The manifolds discussed in the text are best visible on 
the very right as quasi-equidistant bunches of states, (c) Nor- 
malized density of avoided crossings p ac (F) in the vicinity of 
F r =i as a histogram and a Gaussian fit to it (dashed line). 



where the bosonic annihilation (creation) operators at 
the Z-th site are ai(aj), and the number operators are 
nf — ajai, with band index a. J a =a,b are the hop- 
ping amplitudes, and the on-site interparticle interac- 
tion per band has a strength W a . The coupling be- 
tween the bands is given by the repulsive interaction 
W x at simultaneously occupied sites (see Fig. pja)), and 
by the dipole-like term with strength FCo jl4| . The 
latter coupling modifies the above resonance formula 
to F r = A g /2i\^/r 2 — 4Cq by taking into acc ount the 
Stark shift of the levels [14]. The prefactors in (Al) are 
computed based on the single-particle Wannier functions 
ipa(%) — (x\^a) localized within each site (i.e. within 
each double well), given the band structure parameters 
Vo and so [11] . Imposing peri odic boundary conditions in 
space, the Hamiltionian ( |A1| ) is translationally invariant. 
In addition, it is periodic in time, H(t + Tg) = H(t), 
with the Bloch period T# = 1/F. This allows us to 
work with the Floquet Hamiltionian Hf = H(t) — idt, 
which is numerically diagonalized by a Lanczos algo- 
rithm [19J. This procedure is equivalent to diagonal- 
ize the evolution operator integrated over one Bloch pe- 
riod U (Tb, 0) = T exp -i J q Tb H(t)dt , where f is the 
time ordering operator. The quasienergies (eigenvalues of 
Hf) lie within the so-called Floquet zone (FZ) of width 
uob = 2i\F . The extended spectrum is thus given by the 
operation e — >• e + nco>#, with the index n of the FZ. 
For TV atoms, the number of quasienergies is given by 




FIG. 3: (Color online) Regular-to-chaotic transition: The 
main panel shows x as a function of the filling factor N/L. 
Inset: P(s) before (red/grey solid line: N/L = 4/7) and after 
(black solid line: N/L = 9/4) the saturation of %. The dashed 
lines show the Poisson (red/grey) and Wigner-Dyson (black) 
predictions respectively. The parameters are Vb = 5 = so and 
F r = 2 0.0124. 



M s = (N + 2L - 1)\/[LN\(2L - 1)!], considering the re- 
duction by a factor L arising from the periodic boundary 
conditions [TP] , yet the effective dimension of Hf is much 
larger: M s Ak : with Ak = 10 ... 50 Floquet blocks. 

The force is our control parameter to analyze the spec- 
trum in the plane e-F as shown in Fig. [2|a) for the single 
particle case. The gap A gi typically the largest energy 



scale in Eq. (Al) for a realization as shown in Fig. [I] 
allows us to split up the spectrum into subsets of states, 
each labeled by the number of particles in the upper band 
A/5. In the off-resonant regime, where F is not close to 
F r , these subset manifolds do not mix and the levels of 
the various manifolds move parallel with F. Here we have 
approximately H « ©^=0^/5 with an energetic separa- 
tion between the lowest states in energy from the nearest 

manifolds given by A r - A^(l - F/F r ) 2 + 4F 2 C$ / A 2 g 

(see appendix). These manifolds mix, however, strongly 
in the vicinity of F r , where the spectrum shows a large 
number of avoided crossings, see Fig. [2|c). The density 
of avoided crossings p ac (F) sensitively depends on the 
filling factor of the lattice N/L and the interband re- 
pulsive interaction W x . The effective region of manifold 
mixing, the RET regime, is that range of F for which 
the area under the curve p ac (F) is visibly different from 
zero. There A/5 is no longer a good quantum number as 
shown in Fig. [2^b) . Here even the identification of the 
otherwise most distant manifolds ^5 = and A^ = N is 
difficult due to the strong level mixing. 

Regular-to- chaotic transition. The complexity of the 
spectra is now further investigated by means of the 
level spacing (or local gap) distribution P(si) with Si = 
Ei+i — Ei, where (si) = 1, after an appropriated unfold- 
ing procedure [12J [20] . Squeezing more many-body lev- 
els into a single FZ leads to a crossover between regular 
(Poisson) and quantum chaotic (Wigner-Dyson) statis- 
tics of P(s). The number of levels can be increased by 
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varying the filling factor N/L, and the width of the FZ 
by the order r of the resonance. The parameter so con- 
trols the single particle band gap and the manifold split- 
ting A r . As the manifolds approach each other more and 
more in the RET region, nontrivial level statistics are ob- 
served in the range AE of the spectrum, with the many 
avoided level crossings around F r (see Fig. [2jb)) . Strong 
mixing of levels and long-range spectral correlations oc- 
cur as soon as AE becomes of the order of ujb as N/L 
increases. Here we observe Wigner-Dyson statistics in 
P(s), see the inset in Fig. |3j and a logarithmic scaling 
of the level-level correlation function ^ 2 [20] (not shown 
here). To estimate the critical parameters, at which the 
crossover is reached, we propose to study the distribution 
of the energies in AE by means of the quantity 
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with p(e) being the level density of the rescaled ener- 
gies e' — >• 2s/uuB' X characterizes the distribution of 
the levels in energy space. This function saturates to 
Xmin = 1/2, if N 3 is large enough and the energies are 
uniformly distributed over the entire FZ. This latter is 
a signature of chaotic spectra. Otherwise \ is larger 
than this saturation value and diverges if AE <C cjg. 
Fig. [3] shows that the saturation value is reached around 
(N/L) c ~ 1. In this spectral analysis we do not take 
into account those systems for which the greatest com- 
mon divisor gcd(7V, L) is larger or equal to one, due to 
the existence of a temporal symmetry of the Hamiltonian 
as reported in [10 . In the following we discuss important 
consequences that emerge from the spectral characteris- 
tics in the RET regime. 

Level diffusion. The structure of avoided crossings 
presented above provides a perfect setup to study dy- 
namical processes arising from a cascade of Landau- 
Zener (LZ) events in the many-body spectrum around 
F r . We first focus on the diffusion process triggered by 
the parametric-time evolution of individual eigenstates 
\£i) with F(t) = F + at, and a = AF/AT. AF is the 
effective extension of the RET regime and AT ^> Tg 
is the time of the evolution from the starting tilt Fo to 
the final value Ff. The temporal dependence of F(t) 
in Eq. ([2| breaks the time periodicity of the Hamil- 
tonian. In consequence, we must compute now the 
temporal evolution explicitly, e.g., by using a fourth- 
order Runge-Kutta method. Let |$ m (0)) be one of 
the eigenstates at fixed Fo = F r — AF/2. To quan- 
tify the diffusion in energy space we compute the pro- 
jections d(t) = (ei(F k )\$m(t)), where \ei(F k )) are the 
eigenstates at the instantaneous tilt F(t) = F k . As a 
function of the local energies, the distribution of coeffi- 
cients {|C^(£)| 2 } is just the local density of states (LDOS) 
Pm(e,t) = \C i (t)\ 2 5(e-e i ) HE]. At F , P m (e,0) is 
^-shaped but as the tilt increases with time, it delocalizes 
because of the many LZ processes. This effect depends 
on a and sensitively on the type of statistical distribution 
of the spectrum close to F r . This is shown in Fig. [4] (right 
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FIG. 4: (Color online) Diffusion: (left panel) Averaged £ _1 (£) 
and S s h(t) for a = 1.12 x 10 -3 . The blue dashed lines rep- 
resent the maximal derealization value predicted by random 
matrix theory (see main text), (right panel) Temporal evo- 
lution of eigenstates of type Nb — 2 and Nb — 3 through 
a Poisson-distributed spectrum (N/L = 5/6) and a GOE- 
distributed one (N/L = 7/4), respectively. The remaining 
parameters are r — 2 for Vb = 5 = sq. 



panels) for two generic spectra corresponding to the fill- 
ings N/L = 5/6 and N/L = 7/4, which are Poisson- 
and GOE-distributed respectively. We characterize the 
spreading over the spectrum quantitatively by the inverse 
participation ratio [21] and the Shannon entropy [9]: 



M s -l 
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(3) 



S Bb (t) = (- logio 1 (A/ - ) £ |a(i)| 2 log 10 (|G(i)| 2 )) (4) 



with the average taken over a large set of similar initial 
conditions. For the Poissonian case we consider the en- 
semble of initial states that lie in the manifold N^ = 2 
(N/L = 6/5), whilst for the GOE case those from the 
manifold N^ = 3 (N/L = 7/4). The number of acces- 
sible eigenstates is A/" s (5/6) = 728 and A/* s (7/4) = 858. 
Complete derealization occurs as the coefficients {|Q| 2 } 
fluctuate around the equipartition condition for thermal- 
ization \Ci(t)\ 2 — 1/J\f s . Here the localization measures 
slowly converge to the values log 10 (£ -1 ) ~ \og 10 (3/Af s ) ~ 
-2.41 and S sh w log 10 (0.48A4)/ log 10 (AT s ) w 0.88 pre- 
dicted by random matrix theory [4!. It is thus clearly 
seen that in the course of the evolution, the system un- 
dergoes a dynamical diffusion (see Fig. [4]) through the 
accessible Hilbert space. This diffusive spread is much 
stronger for a GOE spectrum /i ~ 94% than for a Poisso- 
nian \i ~ 40%, where \i is the percentage of derealization 
at final time AT with respect to the maximal possible de- 
localization defined above. 

Hence the temporal evolution of a given initial state 
when sweeping F(t) across the RET regime sensitively 
depends on the spectral properties of the system. Strong 
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FIG. 5: (Color online) State preparation: Evolution of the 
eigenstate |$o(0)) « 1 11111; 00000) with (a) a/5AF = 50 
and (b) a/5AF = 1.9 from P/ = 0.46 to P = 0.38. (c) 
Success probability Pt ar for the target state |4>/(AT)) = 
|00000; 11111), with N/L = 5/5 (solid line) and survival prob- 
ability of the initial states |$o(0)) « 1 11111; 00000) (dashed 
line). (d) Instantaneous manifold occupation P(Nb) = 
\P Nb \<f>(i)}\ 2 . a/SAF = 1.9 in (c) and (d). 

diffusion along the spectrum occurs for quantum chaotic 
(Wigner-Dyson distributed), for which the instantaneous 
LDOS P m (£, t) is uniform in the entire FZ since |Q(t)| 2 w 
1/A/" S , whilst a full spreading over all levels is prohibited 
by regular (Poisson) spectra. In the latter case, P m (e,t) 
is typically described by the Breit-Wigner formula [22] 
for t > as shown in Fig. [4] (upper right panel). 

State preparation. The knowledge of the diffusion 
properties presented above allows us to prepare specific 
target states with high fidelity. We concentrate now on 
the case of regular spectra since for them a simple linear 
sweep of the tilt can be used to control the respective evo- 
lution. This is a consequence of the coupling of smaller 
subsets of levels as compared to GOE-distributed spec- 
tra. An extension to the GOE case seems possible, yet 
with more complicated pulses F(t) due to the many en- 
ergy scales present in the problem (see also conclusions 
below). The following analysis is made around the reso- 
nance at F r= i = 0.4248 for Vb = 5 and sq = 1. 

Outside the RET regime, our Floquet eigenstates can 
be identified since they essentially correspond to specific 
translationally invariant Fock states in real space. For 
the following examples, we start our parametric evolu- 
tion with initial Mott-insulator-like states of the untilted 
problems: |$ (* = 0)) « | 111...; 000...). The target state 
|$/(AT)) = |000...; 111...) is of particular interest for ex- 
periments, since its preparation is difficult without apply- 
ing additional complications (e.g. photon-assisted transi- 
tions [23 j). Depending on the wanted target state we 
can adapt the rate of the linear sweep F oc t across the 
RET regime. For given F r and AE (see above), the av- 
erage mean level spacing is 5 = AE/N S , which deter- 
mines the time scale of the sweep rate. For a > 1 we 



expect to follow the diabatic arm of the single particle 
avoided crossing in Fig. |2|a). In our many-body setting, 
the corresponding final state is then close to a state of 
the Nt> = manifold, i.e., to a state where all particles 
are in the lower band. On the other hand for a ~ JAP, 
we expect an adiabatic crossing, implying a strong final 
occupation of the high lying manifolds = N — 1 and 
N b = N (see Fig. §b)). 



For the preparation of the different target states we 
use the the following protocol: (i) the system is pre- 
pared in the state |$ (0)) = 1 11111, 00000) in the flat 
lattice with N/L = 5/5. (ii) A sudden sweep (a ^> 1) 
from P = to Ff > F r is applied, which leaves the 
system unchanged with high fidelity. (Hi) This latter 
state is evolved back from Ff to Fq < F r by a con- 
trolled sweep F(t) = Ff — at with relatively small a. 
Finally, (iv), at t = AP, we compute the success proba- 
bility of the target state in the upper band defined by 
Ptar = |(00000; lllll|P(AP,0)|^o(0))| 2 . This proba- 
bility is shown in Fig. [5[c) (solid line), as a function 
of a rescaled by SAF. It is seen that Ptar(<^) has a 
plateau around a c « 25 AP, reaching a value P tar « 0.81. 
The dashed line corresponds to the survival probability 
-fsurv(^) of the initial condition. It is particular illumi- 
nating to use the above introduced manifolds to picture 
the evolution of the many-body states during the sweep. 
Fig. |5|d) nicely shows the transformation of the N b = 
manifold to states living essentially in the N b = N man- 
ifold. Imperfections in the preparation arise particularly 
from the N b = N — 1 contribution visible in the manifold 
distribution at the final time. 



Conclusions. We have studied the spectral properties 
of an experimentally realizable paradigm of a complex 
many-body quantum system, for which the spectra show 
a regular to chaotic transition. The spectral characteris- 
tics can be probed by sweeping across the resonant region 
using the Stark force as a control parameter. This allows 
us not only predictions on the diffusion of our many-body 
system in Hilbert space, yet as an important application 
also the preparation of target states with high fidelity. 
The latter may be increased even further by applying 
methods of coherent control for our many-body quan- 
tum system [24]. In the chaotic regime, such an optimal 
pulse control is definitely needed to steer the state dif- 
fusion. Our results open a route to engineer complex 
many-body systems at will, with immediate experimen- 
tal implications [6j [15j [25] . 
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FIG. 6: Noninteracting many-body spectrum around F r = 
0.1934, at A g = 1.609, Vb = 5 and s = 2, c.f. Fig. 2 in the 
paper. 



Appendix A: Manifolds far from resonance and 
derivation of manifold energy distance 

Here we motivate the manifold energy scaling defined 
in the main text. Let us begin with the noninteracting 
version of the time- independent Hamiltonian, Eq. (1) of 
the paper, 



ff = EE ^-{a] +1 a l + h.c)+Flnt 

a=a,b l — l 
L 

+ [^ g n b l +FC (b}ai + h.c.) 



i=i 



(Al) 



which, for F = 0, is easily diagonalized in momen- 
tum space. The eigenenergies are given by £F=o,A/" b = 
N^Ag + Sq L \ where is the total number of particles 



in the upper band = (X^ n ?). The quantum num- 
bers Nb define sets of M(A^) degenerate states, where 
X^jVfa M{N},) = M S L. These sets we call manifolds. As F 
increases, single energy levels from the different manifolds 
(corresponding to different try to cross at specific 
tilts F r , which defines the resonantly enhanced tunneling 
(RET) regime. Note that the transition between neigh- 
boring manifolds arises from the single particle transition 
term FCq. This latter single-particle process, at RET 
conditions, can be studied by reducing the full Floquet 
matrix to the resonantly coupled terms only, given by an 
effective 2x2 Hamiltonian matrix [26 . Therefore the en- 
ergy differences between the manifolds \Nb) and + 1) 
can be estimated by diagonalizing the corresponding 2x2 
matrix 



2x2 



£f, iVb+i FCo 

FCq £F,N b 



(A2) 



The energies SF,N b are proportional to F by the the 
simple geometic relation €F,N b — £o = Nb F Ag (F r — F) (see 
Fig. [6|. Thus, by diagonalizing i^2x2, we see that the 
energy difference A r shown in Fig. [6] is 



= A 



' £ F,N h ) 



■4FCJ_ 



^(1-F/F r ) 2 +4FC 2 /A2. (A3) 



The interparticle interaction now splits up the internal 
manifold levels and a strong mixing occurs at RET con- 
ditions when the levels come closest, see Fig. [6] here and 
Fig. 2 of the paper. 
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